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Radiative shocks are shocks in a gas where the radiative energy and flux coming from 
the very hot post-shock material are non-negligible in the shock's total energy budget, 
and are often large enough to heat the material ahead of the shock. Many simulations of 
radiative shocks, both in the contexts of astrophysics and laboratory experiments, use a 
grey treatment of radiative transfer coupled to the hydrodynamics. However, the opacities 
of the gas show large variations as a function of frequency and this needs to be taken 
into account if one wishes to reproduce the relevant physics. We have performed radi- 
ation hydrodynamics simulations of radiative shocks in Ar using multigroup (frequency 
dependent) radiative transfer with the heracles code. The opacities were taken from the 
odalisc database. We show the influence of the number of frequency groups used on the 
dynamics and morphologies of subcritical and supercritical radiative shocks in Ar gas, 
and in particular on the extent of the radiative precursor. We find that simulations with 
even a low number of groups show significant differences compared to single-group (grey) 
simulations, and that in order to correctly model such shocks, a minimum number of 
groups is required. Results appear to eventually converge as the number of groups in- 
creases above 50. We were also able to resolve in our simulations of supercritical shocks 
the adaptation zones which connect the cooling layer to the final post-shock state and 
the precursor. Inside these adaptation zones, we find that the radiative flux just ahead of 
the shock in one or several high-opacity groups can heat the gas to a temperature higher 
than the post-shock temperature. Through the use of Hugoniot curves, we have checked 
the consistency of our radiation hydrodynamics scheme by showing that conservation of 
mass, momentum and energy (including radiative flux) holds throughout the computa- 
tional domain for all our simulations. We conclude that a minimum number of frequency 
groups are required to correctly simulate radiating flows in gases whose opacities present 
large variations as a function of frequency. 



1 Introduction 

Radiative shocks are shocks in a gas where the radiative 
energy and flux coming from the very hot post-shock ma- 
terial are non- negligible in the shock's total energy budget 
(Zel'dovich & Raizer 1967; Mihalas & Mihalas 1984). A radia- 
tive precursor is formed ahead of the shock when the forward 
flux of ionizing photons exceeds the flux of atoms approach- 

* Corresponding author. E-mail address: nell.vaytet@ens-lyon.fr 



ing the shock front. These conditions are met when the 
shock velocity exceeds the threshold required to produce the 
necessary photon flux (Keiter et al. 2002). Two regimes of ra- 
diative shocks are often described in the literature. The first 
is called the subcritical regime, where the shock has only 
a transmissive precursor and the temperature just ahead 
of the discontinuity is not equal to the final downstream 
state temperature. The second is known as the supercritical 
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regime which arises as the strength of the shock increases; a 
diffusive region in the precursor appears and the pre-shock 
temperature reaches the final state temperature (see Drake 
2006, for more details) . 

The classical structure of a subcritical radiative shock 
is depicted in Fig. la. The preshock gas is heated by the 
radiative precursor to a temperature T_ and the shock com- 
pression heats it further to a temperature T + which is higher 
than the final post-shock equilibrium state Ti . The gas then 
cools down inside the cooling layer to reach the equilibrium 
state Ti by radiating the excess energy away. The radiation 
is decoupled from the gas inside the cooling layer and the 
transmissive precursor. The pressure gradient in the pre- 
cursor, through the conservation of mass and momentum, 
causes the velocity to decrease and the density to increase 
to a value p_ ahead of the discontinuity. The sharp density 
jump of the shock from p_ to p + then takes place on the gas 
viscous scale. The density increases further from p + to pi in- 
side the cooling layer as the gas contracts (see also Zel'dovich 
& Raizer 1967; Drake 2006). 

In the case of an optically thick supercritical radiative 
shock (shown in Fig. lb), we have T_ « T\. The radiative tem- 
perature is equal to the gas temperature for the most part, 
except that it remains constant across the cooling layer and 
is higher than the gas temperature inside the transmissive 
precursor (Mihalas & Mihalas 1984). As stated above, the 
pressure gradient at the head of the precursor causes the 
density to increase. Since the gas temperature is close to 
being constant in the diffusive part of the precursor, there is 
no more pressure gradient and the density reaches a plateau 
value p_ ahead of the discontinuity. 

The study of radiative shocks begun in the late 1940s 
with theoretical studies on the Rankine-Hugoniot jump con- 
ditions including a non-negligible radiation pressure for very 
energetic shocks (see Blinnikov & Tolstov 2011, for a short 
review). The studies were very rapidly pursued and extended 
in the field of astrophysics since no processes on Earth 
could achieve high enough energies to produce such shocks. 
Radiative shocks are indeed found in novae (Vaytet et al. 
2007; Orlando et al. 2009; Bode & Evans 2008), supernovae 
(Draine & McKee 1993; Ghavamian et al. 2000; Nymark et al. 
2006), stellar atmospheres (Fadeyev & Gillet 2000; Gillet 
2006), accretion processes in star formation (Stahler et al. 
1980; Commercon et al. 2011), symbiotic stars (Falize et al. 
2009; Imamura 1985) and jets (Raga et al. 1999; Reipurth 
& Raga 1999). This omnipresence makes them a key phys- 
ical process at the heart of high energy astrophysics, and it 
is thus essential to fully understand the details of such a 
mechanism. 

In recent years, with the modern advances in technology, 
radiative shocks have been produced in a number of labo- 
ratory laser facilities (see Bosier et al. 1986; Edwards et al. 
2001; Fleury et al. 2002; Reighard et al. 2006; Busquet et al. 
2007; Michaut et al. 2009, for example), where very high- 
energy lasers are used to drive radiative shocks inside gas 
chambers. This allows new diagnostics of the properties of 
radiative shocks with a much more detailed view than 





Figure 1 : Classical structure of a subcritical (a) and super- 
critical (b) radiative shock. The direction of the gas flow is 
from right to left in the frame where the shock is at rest. 
The panels show the gas (solid) and radiative (dotted) tem- 
perature (top) and the gas density (bottom) as a function of 
distance in each case. The position of the temperature and 
density jumps is marked by the vertical dashed line. The 
relative sizes of the layers are for illustration purposes only. 

would ever be possible in astrophysics. Radiative shock ex- 
periments allow for the validation of numerical simulations 
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which are overwhelmingly used in high-energy physics and 
astrophysics to make predictions on high-energy flows. 

Semi-analytic studies of the structure of radiative shocks 
have been carried out by Lowrie & Rauenzahn (2007); Lowrie 
& Edwards (2008), and several comparisons between exper- 
iments and simulations have also been undertaken by Bou- 
quet et al. (2004); Leibrandt et al. (2005); Reighard et al. 
(2007); Gonzalez et al. (2009), for example. One key piece 
of data that is required by the numerical simulations in or- 
der to accurately model the flow is the opacities of the gas in 
which the shock is launched. Gas opacities show large varia- 
tions as a function of temperature and density as well as fre- 
quency, and including detailed opacities in simulations have 
crucial effects on the structures of radiative shocks. Vaytet 
et al. (20 1 1) performed simulations of a radiative shock in Xe 
using a realistic opacity set and showed the importance of 
taking into account the frequency dependence of the opaci- 
ties, rather than simply integrating over all frequencies, as 
is commonly done in simulations of radiative shocks 

This paper aims to build on the idea that a frequency 
dependent treatment of radiative transfer is crucial in sim- 
ulations of radiative shocks. In particular, we experiment 
further with the multigroup method of Vaytet et al. (2011) 
by studying the effect of the number of frequency groups on 
the structures of the radiative shocks (mainly the variations 
in size of the precursor) . We performed simulations of sta- 
tionary radiative shocks (both sub- and supercritical) using 

1 to 100 frequency groups, and the differences between the 
results are discussed. The opacities are a crucial part of the 
radiative transfer model; they govern the amount of energy 
that will be absorbed and emitted by the gas and can hence 
determine the structure of the flow. Laboratory experiments 
use high atomic number gases to launch radiative shocks to 
take advantage of the strong gas heating due to the lower 
heat capacity. Common choices are argon (Ar) and xenon 
inert gases, and in this work we have chosen to use Ar (see 
section 2.5 for more details). 

2 The multigroup RHD simulations of radiative 
shocks 

2. 1 Radiative transfer 

We use the Mi moment model (Levermore 1984; Dubroca 
& Feugeas 1999) to approximate the equation of radiative 
transfer. The Mi method uses the first two moment equa- 
tions governing the evolution of the radiative energy and flux 



dfF 
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where c is the speed of light, a the absorption/emission coef- 
ficient and B the black body specific intensity. E, F, and P are 
the zeroth, first and second moments of the radiation specific 
intensity, namely the radiative energy density, the radiative 
energy flux, and the radiative pressure, respectively. In or- 
der to close system (1), the radiative pressure is expressed as 
a function of the radiative energy and flux following P = BE 
where D is known as the Eddington tensor. The expression 



for D is obtained by minimizing the radiative entropy which 
yields 

1 - v 3v - 1 F«F 

(2) 
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where 
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(3) 
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and / = ||F||/cE is known as the reduced flux. Note that 
by definition of E and F, we have f < 1, which implies that 
the radiative energy is transported at most at the speed of 
light. In one dimension we simply have P — x&- This clo- 
sure relation recovers the two asymptotic regimes of radiative 
transfer. In the free-streaming limit (i.e. transparent media) , 
we have / = 1 and % = 1. On the other hand, in the diffusion 
limit, f — and x — 1/3, which corresponds to an isotropic 
radiation pressure. 

2. 2 The equations of multigroup radiation hydrodynamics 

We use the multigroup version of the Mi model coupled to 
the gas hydrodynamics described in Vaytet et al. (20 1 1) to ac- 
count for frequency dependence of the absorption and emis- 
sion coefficients (see Shestakov & Offner 2008; van der Hoist 
et al. 2011; Zhang et al. 2013, for other examples of Go- 
dunov multigroup methods). In the multigroup model, the 
equations of radiative transfer are integrated into a finite 
number of frequency bins (or groups) and the opacities are 
averaged over the same frequency ranges. The closure re- 
lation 3 is applied within each frequency group. The more 
the frequency groups, the more accurate the methods be- 
comes, but the higher the computational cost. The system 
of multigroup RHD equations in the comoving frame is 
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where c is the speed of light, and p, u, p and e are the 
gas density, velocity, pressure and total energy, respectively. 
Q v is the third moment of the radiation specific intensity; 
the radiative heat flux. Subscripts v denote monochromatic 
quantities, and we also define 

Xg= | X v dv (9) 



which represents for X — E, F, P the radiative energy, flux 
and pressure inside each group g which holds frequencies 
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between fg-1/2 and v g+ i/ 2 . N g is the total number of groups 
and & g (T) is the energy of the photons having a Planck distri- 
bution at temperature T inside a given group. The quantities 
Op g , o Eg and o Fg are the means of the absorption/emission 
coefficient o v inside a given group weighted by the Planck 
function, the radiative energy and the radiative flux, respec- 
tively. The radiative quantities are expressed in the frame 
comoving with the fluid, which allows simple expressions to 
be used for the source terms on the right hand side of equa- 
tions (5) -(8). The terms involving the frequency derivatives 
of the radiative pressure and heat flux d v (v¥ v ) and d v (vQ v ) 
are solved using a finite volume method in the frequency 
dimension (see Vaytet et al. 2011, for details). 

2.3 Numerical method 

We have implemented the multigroup radiative transfer mod- 
ule of Vaytet et al. (201 1) in the 3D radiation hydrodynamics 
second order Godunov code heracles' (Gonzalez et al. 2007). 
It uses an explicit solver for the hydrodynamics and an im- 
plicit solver for the radiative transfer. The Ar gas equation 
of state is a simple modified ideal gas equation of state; the 
ionization energy is neglected but the ionization state is used 
to compute the mean molecular weight which in turn affects 
the gas temperature. The disregard of the ionization energy 
may greately over-estimate the temperature (ionization can 
represent half of the internal energy for mid- to highly ionized 
flows; Drake 2006), but since we focus solely on the differ- 
ences between mono- and multi-frequency radiative transfer 
methods and no comparison between simulations and exper- 
iments is made throughout, this oversight does not matter 
for the purposes of the present paper. 

In the RHD equations, it is not trivial to compute the ra- 
diative energy and flux-weighted mean opacities o Eg and o Fg . 
Common practise is to set o Fg — o Pg and o Fg — o Rg where 
a R is the Rosseland mean opacity. In this work, we have 
used an average opacity o Eg — o Fg — o Ag which varies be- 
tween the values of o Pg and o Rg depending on the reduced 
flux / (see Appendix A for more details) . However, we wish 
to point out that the inaccuracies which arise from these 
different approximations are reduced as the number of fre- 
quency groups used increases, since in the limit of infinite 
frequency resolution, all of these quantities simply reduce to 
o v . Any approximation is thus less crude in a multigroup 
model than in a grey model. The simulations were run on a 
varying number of CPUs ranging from 12 for the low num- 
bers of frequency groups to 200 for the heavier calculations. 

2. 4 Initial and boundary conditions 

The simulations of stationary radiative shocks were per- 
formed in a one-dimensional regular cartesian grid compris- 
ing 1000 cells (see Appendix B for a discussion on resolu- 
tion) . The grid sizes were L = 1 cm and L = 6 cm for the sim- 
ulations of subcritical and supercritical radiative shocks, re- 
spectively. The discontinuity was initially located at x s — L/4 

^http://irfu.cea.fr/Projets/ Site_heracles / 
*http://irfu.cea.fr/Projets/ Odalisc/ 



Table 1 : Simulation results 



Run 


Number 


Shock 


Mach 


Shock 


Precursor 


name 


of 


velocity 


number 


position 


size 




groups 


(km s _1 ) 




(cm) 


(cm) 


SUB001 


1 






0.250 


0.016 


SUB005 


5 






0.250 


0.025 


SUB010 


10 


30 


5.62 


0.250 


0.024 


SUB020 


20 


0.250 


0.052 


SUB050 


50 






0.250 


0.070 


SUB100 


100 






0.250 


0.073 


SUP001 


1 






1.269 


2.277 


SUP005 


5 






1.245 


2.475 


SUP010 
SUP020 


10 
20 


100 


18.75 


1.239 
1.233 


2.582 
2.650 


SUP050 


50 






1.227 


2.790 


SUP 100 


100 






1.227 


2.811 



Note: The position of the shock is defined as the position where the deriva- 
tive of the velocity is the maximum. The size of the precursor is measured 
between the shock position and the first point (from the right hand side) 
where the gas temperature exceeds 1 . 1 eV. 

and the gas to the right of the discontinuity was given a 
density of p = 10~ 3 g crrT 3 and a temperature of k B T = 1 
eV. The radiative temperature was in equilibrium with the 
gas and the radiative flux was zero. The upstream veloc- 
ity was set to Uo = -30 km s _1 in the subcritical case and 
ito = -100 km s 1 for the supercritical shock. Once the 
upstream state was chosen, the downstream state was cal- 
culated using the Rankine-Hugoniot jump conditions for a 
radiating fluid (Mihalas & Mihalas 1984), which describe 
the conservation of mass, momentum and energy across the 
discontinuity. We find (using the appropriate states of ioni- 
sation) the downstream quantities for the subcritical shock 
to be pi = 3.65 X 10~ 3 gem -3 , Ui = -8.21 km s 1 and 
k B Ti — 6.92 eV. In the case of the supercritical shock, we 
obtain pi = 3.97 X 10~ 3 gem 3 , Ui = -25.18 km s 1 and 
k B T\ — 32.85 eV. The upstream and downstream values are 
also imposed inside ghost cells at the right and left bound- 
aries of the computational domain, respectively. In comput- 
ing the upstream and downstream states, we have assumed 
that we are sufficiently far from the shock so that the radia- 
tive temperature is in equilibrium with the gas temperature 
and that the radiative flux is zero, which is the case in our 
simulations. 

2. 5 The Argon opacities and the decomposition of the fre- 
quency domain 

The opacities for the Ar gas were taken from the odalisc 1 
database, which provides spectral opacities as well as mean 
opacities (Rosseland and Planck) of many elements for a wide 
range of physical conditions. We used the Ar opacities in 
the frequency range hv = 1 — 16,000 eV, computed with 
the potrec code (Mirone et al. 1997) which is based on the 
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Figure 2: The Ar opacities with the decomposition of the frequency domain for 1 to 100 groups. The black line represents 
the opacities of the gas in the pre-shock (p, T) state while the magenta line is for the post-shock state of the supercritical 
case. The colour-bar at the bottom of each frame codes for the group number; this is used in Figs 3, 5 and 7. 



average atom model, including ^splitting and An — transi- 
tions. The spectral opacities for two different densities and 
temperature are shown in Fig. 2; the pre-shock state is rep- 
resented by the black line (p = 10~ 3 g cm" 3 ; k B T = 1 eV) 
while the magenta line is for the supercritical post-shock 
state (p, ~ 4 X 10~ 3 g cm" 3 ; k B Ti ~ 33 eV). One can see that 
the opacities show important orders-of-magnitude variations 
as a function of frequency as well as gas density and tem- 
perature, and this can have a large impact on the results. 
A frequency-averaged model is in this case not an accurate 
approximation since it cannot model situations where a gas 
would be optically thick in one part of the spectrum and 
optically thin in another. 

In this work, we performed simulations with 1,5, 10, 20, 
50 and 100 frequency groups. Choosing the appropriate de- 
composition of the frequency domain among the groups is 



not very straightforward. Ideally, as soon as a large varia- 
tion in k v as a function of frequency is present, one would 
require a new frequency group. When only a small num- 
ber of groups is used, different boundary choices can have 
different impacts on the simulation results. Group bound- 
ary placing naturally becomes less and less important as the 
number of frequency groups increases. For the benefit of a 
fair comparison between simulation results, the decomposi- 
tion of the frequency domain among the groups was simply 
done logarithmically, as shown in Fig. 2. 

3 Results 

We performed simulations of subcritical and supercritical ra- 
diative shocks, using 1 - 100 frequency groups in both cases. 
The properties of the different runs are listed in Table 1. 
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Figure 3: Stationary subcritical radiative shocks using (from top to bottom) 1 to 100 frequency groups. From left to right: 
gas density, gas (black) and radiative (colours) temperatures, radiative flux and gas opacity. At the end of each row is a 
colour bar coding for the different frequency groups. In the two central columns, the magenta curves represent the sum over 
all groups for the radiative temperature and flux. Dashed lines represent negative values of the radiative flux (i.e. flowing 
from right to left) . 



The initial discontinuity, set up with the Rankine-Hugoniot 
conditions described in 2.4, was left to evolve until all the 
structure in the radiative shock was fully developed and the 
stationary regime was reached. All the results shown be- 
low (apart from figures showing a time evolution) are in the 



stationary regime. 

3. 1 The subcritical case 

The results for the simulations of a subcritical radiative 
shock using 1-100 frequency groups are shown in Fig. 3. 
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Figure 4: Left column: gas temperature as a function of distance using 1-100 frequency groups (see color legend at the top) 
in the case of subcritical (b) and supercritical (d) radiative shocks. Right column: distance between the head of the radiative 
precursor and the shock density jump in the case of subcritical (a) and supercritical (c) radiative shocks. 



Each row is for a different number of groups. The columns 
from left to right display (as a function of distance) the gas 
density, the gas (black) and radiative (colours) temperatures, 
the radiative flux and the gas opacity. In the two middle 
columns, the magenta curves represent the sums over all 
groups. 

We first take a look at the mono-group simulation; run 
SUB001 (top row). The profiles exhibit the classic structure 
of a subcritical radiative shock; a transmissive radiative pre- 
cursor extends ahead of the shock heating the upstream gas 
and altering its density and velocity (not shown), the radia- 
tive flux (c) is maximum at the density jump (position of the 
hydrodynamic discontinuity) and we also note the presence 
of a cooling region (or Zel'dovich spike) in the temperature 
plot (b) at the position of the density jump, where the gas 
temperature exceeds the post-shock temperature. The pre- 
cursor measures approximately 0.02 cm and the gas opacity 



(which is averaged over the entire frequency range) lies be- 
tween 5 X 10 4 and 3 X 10 5 cm 2 g _1 throughout (note that the 
size of the precursor is measured between the shock posi- 
tion and the first point from the right hand side where the 
gas temperature exceeds 1.1 eV). 

The first multigroup simulation was performed using 5 
frequency groups (SUB 00 5), the results of which are shown 
in the second row. The first major difference that emerges 
when compared to the grey simulation is that the extent of 
the radiative precursor has increased; even though the gas 
temperature profile has not changed significantly, the radia- 
tive temperature curve extends much further ahead of the 
discontinuity. This can be explained by the range of opac- 
ities observed in the different groups (see Fig. 3h), and in 
particular the opacity of group 3 (light green), whose radia- 
tive energy dominates inside the precursor, which is almost 
an order of magnitude lower than the grey average. The grey 
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Figure 5: Same as Fig. 3 for the supercritical shock case. 



opacity is biased towards a higher value which greatly affects 
the results, especially the radiative flux which now shows a 
long tail extending towards the right end of the simulation 
box. 

The results of the subsequent runs, using 10, 20, 50 
and 100 groups, are shown in the lower rows. While runs 
SUB005 and SUB010 are very similar to each other, the gas 
temperature profile has noticeably changed further in runs 
SUB020 to SUB100, this being most visible between 0.25 



and 0.3 cm where the gas temperature is significantly al- 
tered by the radiation. The gas temperatures are overlayed 
inside the same window in Fig. 4a for a clearer view. The 
position of the density jump (or hydrodynamic shock) and 
the size of the precursor for all the simulations are listed 
in Table 1. The precursor sizes keep increasing with the 
number of frequency groups used. A greater resolution in 
the frequency domain allows an accurate treatment of the 
sharp slopes in the opacity curve (cf. Fig. 2), which affects 
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the amount of absorption ahead of the shock and defines the 
extent of the precusor. The extent of the precursor is plotted 
as a function of time in Fig 4b. The influence of the number 
of groups described above is very clear on this plot. 

3.2 The supercritical case 

The results for the simulations of a supercritical radiative 
shock using 1 - 100 frequency groups are shown in Fig. 5. 
The different rows and columns are identical to Fig. 3. The 
profiles are this time characteristic of a supercritical radia- 
tive shock; a diffusive radiative precursor strongly heats the 
gas ahead of the shock, notably altering the gas density, 
the post-shock gas and radiative temperatures are identi- 
cal to the pre-shock ones and the Zel'dovich spike is clearly 
visible. As in the subcritical case, the size of the precusor 
increases with the number of frequency groups, growing by 
more than 20% between 1 and 100 groups, as explicited in 
Table 1 and illustrated in Figs 4c and 4d. The fact that there 
is very little difference between the 50 and 100-group simu- 
lations indicate that we have almost reached convergence of 
our results. The number of frequency groups used impacts 
the results less than in the subcritical case, since most of 
the radiative precursor is in the diffusive regime where the 
grey approximation is deemed to be accurate. 

These findings have important consequences on predic- 
tions made by numerical simulations on the structure of 
radiative shocks. Indeed, due to the very different precursor 
sizes, studies of radiative shocks which make use of com- 
parisons between observations and numerical calculations 
will most probably be inaccurate if a grey radiative transfer 
model is used. One can compare the right column of Fig. 4 
to the shock-precursor position diagrams found in the liter- 
ature (see Fig. 3 in Michaut et al. 2009 and Figs. 6 and 7 in 
Gonzalez et al. 2006 for example). 

3.3 Electron densities 

In order to illustrate the differences between the grey and 
multigroup simulations further, we now consider the effects 
of our results on the observables commonly obtained in ra- 
diative shock experiments. We compare in Fig. 6 the evo- 
lution of the electron density JV e as a function of time and 
distance in the grey and the 100-group simulations for the 
subcritical (left column) and the supercritical shocks (right 
column). The middle panels (b), (c), (g) and (h) show the N e 
distribution in the simulation frame. The dark red region is 
the post-shock final state with a dense and highly ionised 
medium, while the dark blue region is the pre-shock initial 
state. Between the two, the shock precursor is clearly visible 
in yellow. The sharp transition between the precursor (yel- 
low) and the post-shock state (red) represents the position 
of the density jump as a function of time. The simulation 
frame panels are very useful in highlighting the differences 
between the grey and multigroup simulations; the radiative 
precursor is unmistakably larger in the multigroup case. The 
top panels (a) and (f) show a slice extracted from the simu- 
lation frame data, for an alternative view. We note that for 
the supercritical shock, the rightmost tip of the precursor is 



much sharper in the grey than in the multigroup case. The 
supercritical simulation frame figures also reveal the slight 
displacement of the density jump towards the left, as the 
whole structure of the radiative shock develops over time. 

The bottom panels (d), (e), (i) and (j) show the N e distri- 
bution in the laboratory frame, and are reminiscent of Fig. 3 
in Michaut et al. (2009). These simulated laboratory-frame 
diagnostics show that the differences between the grey and 
multigroup simulations would be large enough to be detected 
in experimental observations. 

3. 4 Detection of adaptation zones around the Zel'dovich tem- 
perature spike 

We now turn to describe what is happening in the vicinity 
of the hydrodynamic discontinuity. Figure 7 shows a close- 
up around the Zel'dovich spike for runs (from left to right) 
SUP 001, SUP 005 and SUP 100. The temperature profile of 
run SUP001 shows that the classic cooling layer structure, 
with a sharp edge on the right and a smooth cooling (or re- 
laxation) region on the left (as depicted in Fig. ??), is almost 
resolved (see Appendix B for a discussion on resolution). 

However, the temperature profile of run SUP005 shows a 
rather different structure, especially on the right hand side. 
There is a smooth region of 'over- temperature' to the right 
of the spike, spanning from x — 1.26 to 1.34 cm (indicated 
by arrow Al), where the gas temperature T_ is above the 
radiation temperature and higher than the final equilibrium 
temperature T\ . It progressively returns to the same value as 
the radiation temperature as we move away from the spike 
towards the right. In a supercritical radiative shock, the pre- 
shock gas temperature cannot in principle exceed the post- 
shock temperature (Zel'dovich & Raizer 1967). The solutions 
are converged in both cases and the radiative precursors do 
not extend out to the right edge of the box. This effect must 
have another origin. 

Looking at the radiative fluxes (panels g and h), we can 
see that in run SUP 1 , the flux shows a plateau right after 
the discontinuity. In SUP 00 5, the same is observed for the 
dominant group 3 (light green) but group 4 (orange) does not 
show a plateau to the right of the shock, instead it decreases 
smoothly as we move away from the shock (see arrow A2). 
This indicates that the radiation in group 4 is being absorbed 
by the material to the right of the shock. Indeed, the group 4 
opacity of the pre-shock material (panel n) is about an order 
of magnitude higher than the grey average in SUP 1 (panel 
m) where the radiation can escape 'freely'. 

It appears that the high opacity in group 4 causes the pre- 
shock material to absorb the radiation in that group only, 
and this consequently heats up the gas a little more. It is 
like having a small precursor inside a precursor. This can be 
understood in terms of the following. A small amount of in- 
cident radiation energy absorbed by gas per unit angle, time, 
area and frequency is equal to the specific intensity loss, i.e. 

dEf s = -dl v cos ddAdadt dv (10) 

which for a constant opacity along direction s becomes 

dEf s = k v I v cos 8 dAdadtdv ds . (11) 
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Figure 6: Time evolution of the electron density (JV e ) as a function of distance for the subcritical case (left column) and the 
supercritical case (right column) using 1 and 100 frequency groups. Panels (a) and (f) show a slice through the 1 (black) 
and 100 (red) groups simulation rest frame data. Panels (b) and (g) show the distribution of N e in the simulation frame for 
1 group while panels (c) and (h) are for the 100-group case. The bottom row displays the N e in the laboratory frame. The 
horizontal white line in panels (b), (c), (g) and (h) show the position at which the slices were extracted. The bottom panels 
(d), (e), (i) and (j) show the JV e distribution in the laboratory frame. 
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Figure 7: A close-up view on the region around the Zel'dovich temperature spike for runs SUP001 (left), SUP005 (center) and 
SUP100 (right). From top to bottom: gas density, gas (black) and radiative (magenta) temperature, radiative flux, reduced 
flux and gas opacity. In the temperature panels, the circles indicate the grid cells. As in the previous figures, the magenta 
curves represent the radiative quantities summed over all groups and the dashed lines symbolise negative values of the 
radiative flux. 
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In one dimension, integrating over frequency and solid angle 
this reduces to 



k v J pdydu dv dAdtds 

/"•CO 

= I K V F V 

Jo 



dv dA dt ds 



(12) 



where u is the direction cosine. So the energy absorbed 
is proportional to J k v F v dv. We now suppose that we 
have two frequency groups, and to mimic the situation in 
run SUP 00 5, we define the group quantities in the follow- 
ing way. The two groups have the same width in the fre- 
quency dimension Av, and we further assume the radiative 
flux and gas opacity to be constant within each group, i.e. 

J-oo 
K V F V dv = 2jg KgFgAv. The first group can be compared 

to group 3 in SUP 5; it dominates the radiative energy 

and flux, but has small reduced flux and opacity. On the 

other hand, the second group has a small energy and flux, 

but a large opacity and reduced flux (similar to group 4 in 

SUP 00 5). We choose 
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The integral over the frequency range in equation (12) now 
becomes a sum over the two groups and we have 



E abs oc [jqFj + k 2 F 2 ]Av 

oc [5 x 10~ 3 cE + 5 x 10~ 3 cE ]Ai; 



(13) 



which is twice as much as what is found by considering only 
the dominant group 1 . This shows that a group with very lit- 
tle radiative energy and radiative flux can contribute signif- 
icantly to the total energy absorbed if it has a large reduced 
flux and opacity. The fact that E1/E2 = 100 means that 
including group 2 does not change the total radiative tem- 
perature, but does change the amount of energy absorbed by 
the gas, and the gas and radiative temperatures can hence 
be decoupled. 

The morphology of the Zel'dovich spike changes again for 
run SUP100 (Fig. 7f). The decoupled regions both to the left 
and the right of the hydrodynamic discontinuity are wider 
than for SUP0 05 (see arrows A3 and A4). The radiative flux 
diagram (i) this time shows more components with a notice- 
able downward trend and large reduced fluxes. The range of 
opacities in the numerous groups mean that the radiation in 
different groups are absorbed at different rates, which yields 
a wider relaxation region to the right of the shock. The same 
effect is also observed to the left of the spike where the gas 
and radiation temperatures are also decoupled (with nega- 
tive fluxes travelling from right to left) . 

The structure observed here is what is described as an 
adaptation zone by Drake (2007a,b); a region across which 
the influence of radiation from the cooling layer on the shock 
structure fades away and where the temperature and other 



gas quantities make their final small adjustments in order 
to reach their final steady-state values. Figure 1 in Drake 
(2007a) actually depicts exactly the situation we observe 
here. There are two adaptation zones on each side of the 
cooling layer (i.e. the Zel'dovich spike) where the gas tem- 
perature is higher than the final state Ti downstream of the 
density discontinuity and higher than the precursor tem- 
perature T_ upstream. This is precisely what our simula- 
tion results show for 10 frequency groups and above. Drake 
(2007a) actually depicts the pre-shock temperature just be- 
fore the discontinuity as below or equal to the final state T\ , 
but he does mention that "ongoing numerical work by John 
Castor suggests that the temperature inside the adaptation 
zone, at the actual density jump, may be pulled up above 
[Ti]". 

The extent of the region inside which the radiative flux 
from the cooling layer still has an effect on the surrounding 
gas will inevitably depend on the opacity of the gas, which 
explains the different observed sizes for the adaptation zones 
as we vary the number of groups. The density is also pic- 
tured in Drake's papers as being slightly lower than the final 
state pi inside the downstream adaptation zone and higher 
than the leftmost precursor value p_ in the upstream zone. 
This is also the case in our results, as shown by arrow A5 in 
Fig. 7c. 

Drake (2007b) provide an analytical estimate of the width 
of the spike for a given shock strength, and when applied to 
our simulation setup, it predicts that the spike should be 
narrower by a factor ~ 30. McClarren & Drake (2010) also 
mention that the Mi model, even though it describes the to- 
tal energy flows correctly, might not perform satisfactorily 
in the vicinity of the cooling layer, but this has yet to be 
investigated. It would appear from the analytical estimates 
that we might not be resolving the real 'physical' Zel'dovich 
spike, but the discussion in Appendix B shows that we are 
sufficiently resolving the spike (for the purpose of this study) 
resulting from our numerical model. We also demonstrate 
that our observation of adaptation zones being absent from 
grey simulations while being clearly detected in multigroup 
simulations is not resolution dependent. While we have to 
acknowledge that we might not be resolving the true width 
of the spike, as we do not attempt to make any comparisons 
with experiments, we believe that our results are still legiti- 
mate. 

Finally, the presence of these adaptation zones in our 
simulations forces us to revise our depiction of a radiative 
shock structure and adopt a more up-to-date description 
(see Fig. 8). The preshock gas is heated by the radiative 
precursor to a temperature T_ ss T\, and increases slightly 
to a temperature T r as we cross the right adaptation zone 
through radiative heating. The shock compression at the 
density jump (discontinuity) heats it further to a tempera- 
ture T + which is higher than the final post-shock equilibrium 
state Ti . The gas then cools down inside the cooling layer to 
reach the intermediate post- shock state T; by radiating the 
excess energy away. The final small adjustments are made 
across the left adaptation zone to reach the final post-shock 
state Ti (see also Drake 2007b). The radiative temperature 
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(dotted line) is equal to the gas temperature for the most part, 
except that it remains constant across the cooling layer and 
the adaptation zones (decoupled from the gas), and is higher 
than the gas temperature inside the transmissive precursor. 
The pressure gradient at the head of the precursor, through 
the conservation of mass (4) and momentum (5), causes the 
velocity to decrease and the density to increase to a value 
p_ ahead of the discontinuity. Since the gas temperature 
is close to being constant in the diffusive part of the pre- 
cursor, there is no more pressure gradient and the density 
reaches a plateau value p_ ahead of the discontinuity. A 
small compression from p_ to p r occurs as we cross the right 
adaptation zone, then followed by the sharp density jump of 
the shock from p_ to p + which takes place on the gas viscous 
scale. The density then rises rapidly to p; inside the cooling 
layer through strong contraction of the radiating gas. Lastly, 
the density slowly reaches the final state pi across the left 
adaptation zone. 




energy hold (see Zel'dovich & Raizer 1967; Mihalas & Miha- 
las 1984, for example). Figure 9 shows (from top to bottom) 
the gas velocity, pressure and radiative flux as a function of 
inverse compression ratio for every grid point in our subcrit- 
ical (left) and supercritical (right) simulations, as well as the 
analytical solutions (solid black line) which are 



u(q) = -Uori 

pin) = Poi^(i -n) + Po 



2y y + 1 2 \ YPoUo 
■n 7U - i|- - — r(! 



y - 1 y - 1 



y- 1 



(14) 
(15) 

n) + F 

(16) 



where y is the ratio of specific heats and the radiative flux F 
is expressed in the laboratory frame (note that we have con- 
verted our radiative quantities to the laboratory frame using 
equations 91.16 and 91.17 in Mihalas & Mihalas 1984). 

This illustrates well how our numerical scheme preserves 
the important physics of conservation of mass, momentum 
and energy. Any simulation of shocks which has points lying 
away from these analytical curves does not conserve these 
fundamental quantities. The points which do not lie on the 
analytical curves in Fig. 9 are all in the shock transition re- 
gion which is inevitably spread over a minimum number of 
cells for a Godunov method and do not strictly match the 
analytical solution (see Sincell et al. 1999; Drake 2007b, for 
example) . This is mostly visible in the subcritical case where 
the radiative flux follows a constant value transition while 
the gas velocity and pressure follow an approximate Hugo- 
niot curve through the transition to connect the upstream 
and downstream states. Note that only a single cell for each 
simulation does not lie on the analytical prediction, illustrat- 
ing the fact that the transition is spread over only one or two 
grid cells. We take the opportunity here to point out that 
the comoving frame formalism, which has been criticized in 
several works for not rigorously conserving the total energy 
(see Mihalas & Klein 1982; Krumholz et al. 2007, for ex- 
ample), appears here to perform very well in conserving the 
fundamental quantities. 

4 Conclusions 



Figure 8: Structure of a supercritical radiative shock 
(adapted from Drake 2007a). The direction of the gas flow 
is from right to left in the frame where the shock is at rest. 
Top panel: gas (solid) and radiative (dotted) temperature as 
a function of distance. Bottom panel: gas density as a func- 
tion of distance. The position of the temperature and density 
jumps is marked by the vertical dashed line. The relative 
sizes of the layers are for illustration purposes only. 

3.5 Hugoniot curves 

In this section we look at the Hugoniot curves for gas ve- 
locity, pressure and radiative flux in our simulations. The 
Hugoniot curves are analytical predictions for the state of 
gas quantities as a function of the inverse compression ratio 
U — Pal P f° r which conservation of mass, momentum and 



We have performed simulations of stationary radiative 
shocks in Ar gas using a multigroup radiation hydrody- 
namics scheme. Gas opacities depending on temperature, 
density and frequency were used in the equations of ra- 
diation hydrodynamics to achieve convincing results. The 
simulations reproduced all the detailed structure of a ra- 
diative shock, including the radiative precursor, the cooling 
layer and even the adaptation zones connecting the cooling 
layer to the final downstream state and the precursor. Our 
results show that grey simulations produce very different 
results compared to multigroup ones, and that frequency- 
independent calculations are not deemed an accurate de- 
scription of the problem. Indeed, multigroup simulations 
showed increases by a factor of four in precursor size in the 
subcritical case, and an increase of 20% in the supercriti- 
cal case. The simulations with 5 to 100 groups also show 
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Figure 9: Hugoniot curves for the radiative shock simulations using 1 - 100 frequency groups (see legend in the top right cor- 
ner for the meaning of the symbols): gas velocity (top), pressure (middle) and radiative flux (bottom) as a function of inverse 
compression ratio in the case of the subcritical shock (left column) and the supercritical shock (right column) . The analytical 
solutions are overlayed for comparison (black solid line) . The simulation points not lying on the analytical solutions are all 
within the shock transition. 



increasing precursor sizes, with a suspected convergence of 
results for 50 groups and above. Multigroup effects were 
also seen to be important in the vicinity of the cooling layer, 
where adaptation zones absent from grey simulations were 
clearly detected. 

We have to acknowledge that several caveats need to be 
taken into account when considering the results reported in 
this work. Firstly, the exact sizes of the radiative precursors 
are not entirely correct since higher resolution simulations 
reported in Appendix B yield slightly different results (typical 
differences are of the order of 5 - 10%), only the relative in- 
crease in precursor sizes between the different simulations 



are of notorious relevance. Secondly, even though our re- 
sults are numerically quantitatively converged in the prox- 
imity of the Zel'dovich spike, it is not clear whether the M\ 
model performs accurately enough on such small scales. We 
remind the reader that this work is not an attempt at directly 
comparing numerical simulations to experiments, merely a 
study of the effects of frequency dependence on the results 
of radiation hydrodynamic calculations. 

Nevertherless, the findings presented still have important 
consequences on predictions made by numerical simulations 
on the structure of radiative shocks. Indeed, studies of as- 
trophysical (accretion processes, supernova remnants, jets, 
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etc. . . ) or laboratory radiative shocks will most probably be 
inaccurate if a grey radiative transfer model is used. The 
impact can not only be large when looking at the radiative 
precursor sizes, but also on the total energy budget, deter- 
mining the amount of energy converted into radiation and 
absorbed by pre-shock material. 

It is difficult to compare the results of this work with ex- 
perimental data directly since our idealised setup of station- 
ary radiative shocks is very far from the situation in labora- 
tories where laser-driven radiative shocks travel through gas 
chambers and it is often unclear if they ever reach a station- 
ary state. Realistic calculations using a piston-like shock- 
driving boundary, as well as a more sophisticated equation 
of state, will be more appropriate for conducting detailed 
modelling of laboratory experiments. 
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A Using the correct opacity average 

In the RHD equations (5)-(8), it is challenging to compute 
the radiative energy and flux-weighted mean opacities k e 
and k f (note that a — Kp), since the quantities E v and F v 
are not necessarily known at all wavelengths. The choice 
of approximation for these quantities, which are crucial to 
the RHD calculations, is not trivial (see Pomraning 1973, p. 
88+ for a discussion). Common practise is to set k e — k p and 
k f — k r where k p and k r are the Planck and Rosseland mean 
opacities, respectively (see Larsen & Lane 1994; Offner et al. 
2009, for example). However, for the Mi model, numerical 
stability is immensely improved if an identical value for k e 
and Kp is used, which prompts us to use an average value 

Ka V = K E = Kp. 

The Planck and Rosseland means are applicable to dif- 
ferent regimes. In the diffusion limit, the opacity should 
be equal to the Rosseland mean, while the Planck mean is 
appropriate in the free-streaming limit. We need to make 
sure that we recover this behaviour with our average opacity. 
Sampson (1965) proposed an average which varies between 
K P and Kr depending on the optical depth. Here, we use the 
reduced radiative flux / = ||F||/cE as a measure of the dif- 
fusivity (optically thick) or transmissivity (optically thin) of 
the flow. We define a parameter a(f ) which varies between 
zero and one according to f which then allows us to write 
the average opacity 

Ka V = (1 — a)K R + OKp. (17) 

In order to recover the diffusion and free -streaming limits, a 
needs to have the following properties: a — > when / — > 
and a — > 1 when / — > 1 . The simplest formula with these 
properties is just the linear function a — f. However, we 
argue that the a parameter is meant to represent the tran- 
sition from a regime dominated by diffusion to a radiation- 
dominated regime. In this sense, one expects the transition 
from one to the other to be rather rapid, as opposed to a 
smooth linear averaging between the k p and k r values. We 
thus propose a different formula for a which will better re- 
produce this behaviour. After some experimenting, we finally 
chose 



which is plotted in Fig. 10 (red) alongside the simpler a — f 
(black). We would like to point out here that we have no 
physical explanation for equation (18), it is simply an ad-hoc 



choice of a function with the correct properties. However, 
we can justify our choice by testing the method in a simple 
case. 




/ 



Figure 10: a parameter as a function of reduced flux/: the 
black solid line is the simple a—f while the red line is given 
by equation (18). 

We performed eight simulations with a single frequency 
group; four calculations of a subcritical radiative shock and 
four of a supercritical shock. In each case, the first was 
carried out using a—f, the second using a — a s , the third 
using a — (which corresponds to k^ v — k r ) and the fourth 
using a = 1 (Ka V = k p ) . We show the results in Fig. 1 1 . 
In the case of the subcritical radiative shock (left column), 
we see that the average opacities (black and red) are equal 
to the Rosseland mean far away from the discontinuity (on 
both sides). They then become close to the Planck average 
in the intermediate region (between 0.25 and 0.28 cm) where 
the reduced flux is large. The temperature plot shows that 
the results for both the averaging schemes are very close to 
the blue curve (Planck mean) which is the desired result in 
this optically thin regime. For the supercritical shock (right 
column), the average opacity should produce results which 
resemble the green curves, which are appropriate in this dif- 
fusive regime (only a small region of the grid has a large f] . 
This time, the simple a—f approximation shows its limita- 
tions with a precursor size about two thirds of the size of the 
one observed in the Ka V = k r case. On the other hand, the 
a — a s model performs much better, producing results very 
similar to the k^ v — k r simulation. We thus conclude that 
the expression given in equation (18) is an effective model 
to simulate problems in both diffusive and free-streaming 
limits. 

In a final note, we wish to reiterate that the inaccuracies 
which arise from the approximation of setting k e — Kp — 
are reduced as the number of frequency groups used in- 
creases, since in the limit of infinite frequency resolution, 
all of these quantities simply reduce to k v . This approxima- 
tion is thus less crude in a multigroup model than in a grey 
model. 
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Figure 11: Comparison of the different opacity averaging functions. The left and right columns show the results for the 
subcritical and supercritical radiative shocks, respectively. The top row displays the gas (solid) and radiative (dashed) tem- 
peratures, the middle row is the gas opacity and the bottom row shows the reduced flux. The colour-coding is as follows: 
black is for a —f, red is for a given by equation (18), green is for Ka V = k r and blue is for Ka V = k p . 



B A comment on spatial resolution 

In Fig. 7 it would appear that with our spatial resolution 
of 1000 grid cells, we do not fully resolve the very narrow 
Zel'dovich spike. For completeness, we report in this section 
a short spatial resolution study to confirm that our results 
concerning the influence of the number of frequency groups 
on the size of the radiative precusors and adaptation zones 
on each side of the spike are robust. Figure 12 compares the 
effects of grid resolution and number of frequency groups on 
the shock structures. The top row shows the temperature 
profiles in the vicinity of the cooling layer for 1 frequency 
group using 500, 1000, 5000 and 10,000 cells in the sub- 
critical (a) and supercritical (b) cases (see color key at the 
top of the figure). In the subcritical case, the different reso- 



lutions appear to yield similar results. For the supercritical 
runs, the spike appears well resolved in the high resolution 
simulations. It looks thinner than in the lower resolution 
calculations and displays a (smooth) maximum to the left 
of the discontinuity, a structure that much resembles the 
results of Lowrie & Edwards (2008). At a Mach number of 
~ 19, the spike structure fits well between their depictions 
of shocks at M = 3 (Fig. 10) and M = 27 (Fig. 13). A strong 
convergence of results is observed between 5000 and 10,000 
cells. 

The middle row shows again the temperature profiles of 
the cooling layer for different resolutions but this time using 
5 frequency groups. Panel (c) reveals that 500 cells is most 
probably too few to resolve the cooling layer. Both panels 
also show convergence between 5000 and 10,000 cells. More 
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Figure 12: Comparison between the effects of grid resolution and number of frequency groups. Top row: temperature 
profiles in the vicinity of the cooling layer for 1 frequency group using different spatial resolutions in the subcritical (a) 
and supercritical (b) cases (see color key at the top of the figure). Middle row: same as the top row but using 5 frequency 
groups. Bottom row: size of precursor as a function of time for the 1-group (dashed) and 5-groups (solid) simulations in the 
subcritical (e) and supercritical (f) cases. 



importantly, panels (b) and (d) demonstrate that our over- 
all conclusions on the presence of adaptation zones in the 
multigroup simulations remain. Even though the relaxation 
region to the right of the spike is more pronounced in the 
5000-cell simulation (grey), the results from the 1000-cell 
simulations (cyan) agree qualitatively; an adaptation zone 
is absent from the grey simulations but is detected in the 
multigroup simulations. 

Analytical estimates of the width of the spike from Drake 
(2007b) suggest that the real physical spike for the same val- 



ues of shock velocity and initial state density might in fact be 
much narrower (by a factor of ~ 30). Nevertherless, the fact 
remains that we have converged spatially on the structure of 
the cooling layer given by our Mi model, which is what mat- 
ters for the present study. In addition, analytical estimates 
also make use of approximations and the true width of the 
spike is probably not known accurately. 

The bottom row displays the size of the radiative pre- 
cursors as a function of time for the different simulations. 
To better distinguish the separate runs, we have plotted the 
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simulations using 1 group with dashed lines and the simula- 
tions using 5 groups with solid lines (the colours remain the 
same as in the other panels). These plots reveal that the spa- 
tial resolution also affects the total size of the precursor (here 
again the results have converged for 5000 cells and above). 
However, as we do not make any direct comparisons with ex- 
periments or observations throughout this work but are only 
interested in the relative differences in precursor sizes, we 
argue that our conclusions regarding the increase in precur- 
sor size as a function of number of frequency groups remain 
qualitatively correct. Moreover, panels (e) and (f) show quite 
clearly that the difference in precursor size due to a change 
in number of groups is larger than the difference observed 
from a change in number of cells. In addition, Fig. 12 only 
shows the results using 1 and 5 groups, and the size of the 



precursor is seen to continue increasing all the way up to 100 
groups (see Fig. 4). In the case of the subcritical shock, dif- 
ferences between runs with 1000 and 10,000 cells are ~ 20% 
and ~ 8% for the 1 -group and 5-groups simulations, respec- 
tively, while differences between the 1 -group and 100-groups 
(1000 cells) simulations are higher than 400% (see Table 1). 
As for the supercritical, resolution alters the precursor sizes 
by only 2 — 4% while frequency groups have an effect of the 
order of 20%. Finally, the relative differences in precursor 
size between 1- and 5-group simulations for a given reso- 
lution remain approximately constant. It was not possible 
for us to run a 100-group simulation using 5000 cells on 
a realistic timescale, but we believe that the results of the 
1000-cell simulations are qualitatively robust. 



